Progress in the R ecosystem for open source spatial analysis software
Roger Bivand
Thursday 12 December 2019, 16:45-18:45, sala 108
Copyright
All the material presented here, to the extent it is original, is available under CC-BY-SA. Parts build on joint tutorials with Edzer Pebesma.
Required current contributed CRAN packages:
I am running R 3.6.1, with recent update.packages().
needed <- c("raster", "rgdal", "RSQLite", "stars", "abind", "sp", "mapview",
"sf", "osmdata", "wordcloud", "RColorBrewer")Script
Script and data at https://github.com/rsbivand/uw_191212/raw/master/uw_Rgeo_191212.zip. Download to suitable location, unzip and use as basis.
Script with local examples: https://github.com/rsbivand/uw_191212/raw/master/thurs191212.R
Preface/positionality
I am afraid that this talk will not contain many graphical elements. Why is this? Text is processed by different parts of our cognitive systems than images, which open for many more mis-understandings. Text is generally a more direct channel of communication.
Code is privileged text, it shows what is happening explicitly. It facilitates the reproduction of the workflow, and enables clear checking of robustness. Literate programming, and its modern form as notebooks that can be executed to reproduce outcomes, may seem hard to digest, but they enable the reader to become a co-creator of content.
It is a different kind of performance than a story, in that the reader may enter the workflow. In GUI-based stories, it is very possible that neither the author nor the reader can reconstruct the actual workflow yielding the presented outcomes.
Outline
Introduction and background: why break stuff (why not)?
Spatial data input/output and representation: movement from legacy sp to standards-compliant sf representation; coordinate reference systems; developments similar to GeoPandas and Shapely in Python
Opportunities in visualization (tmap, mapview), but also challenges when upstream software libraries evolve (PROJ/GDAL)
R-Spatial: getting to here and now
Background
In the early and mid 1990s, those of us who were teaching courses in spatial analysis beyond the direct application of geographical information systems (GIS) found the paucity of software limiting.
In institutions with funding for site licenses for GIS, it was possible to write or share scripts for Arc/Info (in AML), ArcView (in Avenue), or later in Visual Basic for ArcGIS.
If site licenses and associated dongles used in the field were a problem (including students involved in fieldwork in research projects), there were few alternatives, but opportunities were discussed on mailing lists.
From late 1996, the R programmimg language and environment began to be seen as an alternative for teaching and research involving spatial analysis.
R uses much of the syntax of S, then available commercially as S-Plus, but was and remains free to install, use and extend under the GNU General Public License (GPL).
In addition, it could be installed portably across multiple operating systems, including Windows and Apple MACOS.
At about the same time, the S-Plus SpatialStats module was published (Kaluzny et al. 1998), and a meeting occurred in Leicester to which many of those looking for solutions took part.
Much of the porting of S code to R for spatial statistics was begun by Albrecht Gebhardt as soon as the R package mechanism matured. Since teachers moving courses from S to R needed access to the S libraries previously used, porting was a crucial step.
CRAN listings show tripack (Renka and Gebhardt 2016) and akima (Akima and Gebhardt 2016) - both with non-open source licenses - available from August 1998 ported by Albrecht Gebhardt; ash and sgeostat (Majure and Gebhardt 2016) followed in April 1999.
The spatial package was available as part of MASS (Venables and Ripley 2002), also ported in part by Albrecht Gebhardt.
In the earliest period, CRAN administrators helped practically with porting and publication.
Albrecht and I presented an overview of possibilities of usin R for research and teaching in spatial analysis and statistics in August 1998 (Roger Bivand and Gebhardt 2000).
The S-PLUS version of splancs provided point pattern analysis (Rowlingson and Diggle 1993, 2017).
I had contacted Barry Rowlingson in 1997 but only moved forward with porting as R’s ability to load shared objects advanced.
In September 1998, I wrote to him: “It wasn’t at all difficult to get things running, which I think is a result of your coding, thank you!”
However, I added this speculation: “An issue I have thought about a little is whether at some stage Albrecht and I wouldn’t integrate or harmonize the points and pairs objects in splancs, spatial and sgeostat - they aren’t the same, but for users maybe they ought to appear to be so”.
This concern with class representations for geographical data turned out to be fruitful.
A further step was to link GRASS and R (Roger Bivand 2000), and followed up at several meetings and working closely with Markus Neteler.
The interface has evolved, and its almost current status is presented by (Lovelace, Nowosad, and Muenchow 2019), we return to the current status below.
A consequence of this work was that the CRAN team suggested that I attend a meeting in Vienna in early 2001 to talk about the GRASS GIS interface.
The meeting gave unique insights into the dynamics of R development, and very valuable contacts.
Later the same year Luc Anselin and Serge Rey asked me to take part in a workshop in Santa Barbara, which again led to many fruitful new contacts (Bivand 2006).
Further progress was made in spatial econometrics (Bivand 2002).
2003 Vienna workshop
During the second half of 2002, it seemed relevant to propose a spatial statistics paper session at the next Vienna meeting to be held in March 2003, together with a workshop to discuss classes for spatial data.
I had reached out to Edzer Pebesma as an author of the stand-alone open source program gstat (Pebesma and Wesseling 1998); it turned out that he had just been approached to wrap the program for S-Plus.
He saw the potential of the workshop immediately, and in November 2002 wrote in an email: “I wonder whether I should start writing S classes. I’m afraid I should.”
Virgilio Gómez-Rubio had been developing two spatial packages, RArcInfo (V. Gómez-Rubio and López-Quílez 2005; Gómez-Rubio 2011) and DCluster (V. Gómez-Rubio, Ferrándiz-Ferragud, and Lopez-Quílez 2005; Gómez-Rubio, Ferrándiz-Ferragud, and López-Quílez 2015), and was committed to participating.
Although he could not get to the workshop, Nicholas Lewin-Koh wrote in March 2003 that: “I was looking over all the DSC material, especially the spatial stuff. I did notice, after looking through peoples’ packages that there is a lot of duplication of effort. My suggestion is that we set up a repository for spatial packages similar to the Bioconductor mode, where we have a base spatial package that has S-4 based methods and classes that are efficient and general.”
Straight after the workshop, a collaborative repository for the development of software using SourceForge was established, and the R-sig-geo mailing list (still with over 3,500 subscribers) was created to facilitate interaction.
Beginnings of sp
So the mandate for the development of the sp package emerged in discussions between interested contributors before, during, and especially following the 2003 Vienna workshop.
Coding meetings were organized by Barry Rowlingson in Lancaster in November 2004 and by Virgilio Gómez-Rubio in Valencia in May 2005, at both of which the class definitions and implementations were stress-tested and often changed radically; the package was first published on CRAN in April 2005.
The underlying model adopted was for S4 (new-style) classes to be used, for "Spatial" objects, whether raster or vector, to behave like "data.frame" objects, and for visualization methods to make it easy to show the objects.
Relationships with other packages
From an early point in time, object conversion (known as coercion in S and R) to and from sp classes and classes in for example the spatstat package (Baddeley and Turner 2005; Baddeley, Rubak, and Turner 2015; Baddeley, Turner, and Rubak 2019).
Packages could choose whether they would use sp classes and methods directly, or rather use those classes for functionality that they did not provide themselves through coercion.
Reading and writing ESRI Shapefiles had been possible using the maptools package (Bivand and Lewin-Koh 2019) available from CRAN since August 2003, but rgdal, on CRAN from November 2003, initially only supported raster data read and written using the external GDAL library (Warmerdam 2008).
Further code contributions by Barry Rowlingson for handling projections using the external PROJ.4 library and the vector drivers in the then OGR part of GDAL were folded into rgdal, permitting reading into sp-objects and writing from sp-objects of vector and raster data.
Completing the sp-verse
For vector data it became possible to project coordinates, and in addition to transform them where datum specifications were available.
Until recently, the interfaces to external libraries GDAL and PROJ have been relatively stable, and upstream changes have not led to breaking changes for users of packages using sp classes or rgdal functionalities, although they have involved significant maintenance effort.
The final part of the framework for spatial vector data handling was the addition of the rgeos package interfacing the external GEOS library in 2011, thanks to Colin Rundell’s 2010 Google Summer of Coding project.
The rgeos package provided vector topological predicates and operations typically found in GIS such as intersection; note that by this time, both GDAL and GEOS used the Simple Features vector representation internally.
Applied Spatial Data Analysis with R (ASDAR, Springer: useR! series)
By the publication of ASDAR (Bivand, Pebesma, and Gomez-Rubio 2008), a few packages not written or maintained by the book authors and their nearest collaborators had begun to use sp classes. By the publication of the second edition (Bivand, Pebesma, and Gomez-Rubio 2013), we had seen that the number of packages depending on sp, importing from and suggesting it (in CRAN terminology for levels of dependency) had grown strongly. In late 2014, (de Vries 2014) looked at CRAN package clusters from a page rank graph, and found a clear spatial cluster that we had not expected. This cluster is from earlier this week:
Spatial data
Spatial data typically combine position data in 2D (or 3D), attribute data and metadata related to the position data. Much spatial data could be called map data or GIS data. We collect and handle much more position data since global navigation satellite systems (GNSS) like GPS came on stream 20 years ago, earth observation satellites have been providing data for longer.
When we began, most vector data had been digitised manually, with manual punching of attribute data, or was shared as downloaded files often as ESRI Shapefiles, the format introduced with ArcView. Nowadays, much data is available from online cloud sources, for example Open Street Map, here used to extract the light rail system in Bergen using the osmdata package. Setting up an Area of Interest (AoI), we query the data source for key values with two qualifications (those uploading the data were not consistent). These are then merged to give a single data object:
suppressPackageStartupMessages(library(osmdata))
library(sf)## Linking to GEOS 3.8.0, GDAL 3.0.2, PROJ 6.2.1bbox <- opq(bbox = 'bergen norway')
byb0 <- osmdata_sf(add_osm_feature(bbox, key = 'railway',
value = 'light_rail'))$osm_lines
tram <- osmdata_sf(add_osm_feature(bbox, key = 'railway',
value = 'tram'))$osm_lines
byb1 <- tram[!is.na(tram$name),]
o <- intersect(names(byb0), names(byb1))
byb <- rbind(byb0[,o], byb1[,o])
saveRDS(byb, file="byb.rds")Spatial vector data is based on points, from which other geometries are constructed. Vector data is often also termed object-based spatial data. The light rail tracks are 2D vector data. The points themselves are stored as double precision floating point numbers, typically without recorded measures of accuracy (GNSS provides a measure of accuracy). Here, lines are constructed from points.
It has also become typical over the past few years to show such data objects in a pan/zoom interactive setting, here with the mapview package:
library(mapview)
mapview(byb)Data handling
We can download monthly CSV files of city bike use, and manipulate the input to let us use the stplanr package to aggregate origin-destination data. One destination is in Oslo, some are round trips, but otherwise things are OK.
bike_fls <- list.files("bbs")
trips0 <- NULL
for (fl in bike_fls) trips0 <- rbind(trips0,
read.csv(file.path("bbs", fl), header=TRUE))
trips0 <- trips0[trips0[, 8] < 6 & trips0[, 13] < 6,]
trips <- cbind(trips0[,c(1, 4, 2, 9)], data.frame(count=1))
from <- unique(trips0[,c(4,5,7,8)])
names(from) <- substring(names(from), 7)
to <- unique(trips0[,c(9,10,12,13)])
names(to) <- substring(names(to), 5)# find origins and destinations and merge them
stations0 <- st_as_sf(merge(from, to, all=TRUE),
coords=c("station_longitude", "station_latitude"))
stations <- aggregate(stations0, list(stations0$station_id),
head, n=1)
suppressWarnings(stations <- st_cast(stations, "POINT"))
st_crs(stations) <- 4326# sum trips between origins and destinations
od <- aggregate(trips[,-(1:4)], list(trips$start_station_id,
trips$end_station_id), sum)We can use the stplanr package to generate the station-to-station desire lines between different stations, omitting bikes dropped at the same station where they were picked up:
# drop equal origins and destinations
od <- od[-(which(od[,1] == od[,2])),]
# create desire lines
library(stplanr)
od_lines <- od2line(flow=od, zones=stations, zone_code="Group.1",
origin_code="Group.1", dest_code="Group.2")
saveRDS(od_lines, "od_lines.rds")Origin-destination lines
od_lines <- readRDS("od_lines.rds")
mapview(od_lines, alpha=0.2, lwd=(od_lines$x/max(od_lines$x))*10)We can use CycleStreets to route the volumes onto OSM cycle paths, via an API and API key.
# allocate desirelines to cycle paths using CycleStreet
Sys.setenv(CYCLESTREET="XxXxXxXxXxXxXxXx")
od_routes <- line2route(od_lines, "route_cyclestreet",
plan = "fastest")
saveRDS(od_routes, "od_routes.rds")Routed lines along cycle routes
od_routes <- readRDS("od_routes.rds")
mapview(od_routes, alpha=0.2, lwd=(od_lines$x/max(od_lines$x))*10)